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A new micro-scale refrigeration system composed of arrays of silicon MEMS cooling elements that 
operate on the Stirling cycle has been designed. In this paper, we describe a multiphysics computational 
approach for analyzing the system performance that considers compressible flow and heat transfer with 
a large deformable mesh. The regenerator pressure drop and effectiveness are first explored to determine 
the optimal porosity. A value near 0.9 is found to maximize the coefficient of performance. To overcome 
the computational complexity brought about by the fine pillar structure in the regenerator, a porous 
medium model is used to allow for modeling of a full element. Parametric studies demonstrate the effect 
of the operating frequency on the cooling capacity and the coefficient of performance. 

© 2013 Elsevier Masson SAS. All rights reserved. 


1. Introduction 

Micro-scale coolers have a wide range of potential applica¬ 
tions, such as cooling chip- and board-level electronics, sensors, 
and radio-frequency systems [1]. The coefficient of performance 
of the cooler (COP), is defined as the ratio of heat removal, Q, to the 
input work, W, i.e., COP = (Q/W). The Carnot coefficient of per¬ 
formance, COPo which represents the maximum theoretical effi¬ 
ciency possible between constant temperature hot (Th) and cold 
( Th ) sources, can be found from COP c = T C /(T H - T c ) 

Thermoelectric coolers are a common refrigeration device and 
have been scaled to the micro-domain [2-5]. The efficiency of 
thermoelectric energy conversion is determined by the thermo¬ 
electric material’s figure of merit, ZT, which is defined by (S 2 crT//<) 
where S, a and k are the Seebeck coefficient, electrical conductivity, 
and thermal conductivity [6]. Significant challenges exist to in¬ 
crease ZT beyond unity due to the difficulty of reducing the thermal 
conductivity while maintaining good electrical properties [7]. The 
maximum COP of a thermoelectric cooler is [2] 
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t c (VTTzT-t h /t c ) 

COP = -^- V (1) 

(T H -T c )(yr+ZT+l) 

This maximum COP is plotted as a function of the temperature 
difference (Th - Tc) in Fig. 1 for ZT values of one and two (the latter 
is the maximum value currently available in bulk thermoelectric 
materials [2]). For a typical temperature difference of 25 K for 
cooling in ambient conditions, the maximum COP of a realistic 
thermoelectric cooler is 2.6, which is 24% of the Carnot COP. 

In recent years, novel solid-state cooling technologies, such as 
electrocaloric and magnetocaloric cooling, have attracted interest 
[8,9]. The electrocaloric effect is a phenomenon in which reversible, 
polarization-related temperature and entropy changes appear un¬ 
der the application and removal of an electric field. Potential 
application of this technology is limited, however, by the relatively 
low temperature and entropy changes at ambient conditions for 
most ferroelectric materials [8]. The magnetocaloric effect is similar 
to the electrocaloric effect, with the use of a magnetic field instead 
of an electric field. The high magnetic fields required for magne¬ 
tocaloric cooling are difficult to realize, especially at the micro-scale 

[ 91 - 

Micro-scale devices operating on the Stirling cycle, which was 
developed in 1816 [10], are an attractive potential alternative due to 
the high efficiencies realized for macro-scale Stirling machines [11 ]. 
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Nomenclature 

P 

pressure, Pa 



Q 

cooling capacity, W/cm 2 

COP 

coefficient of performance 

t 

time, s 

COP c 

Carnot COP 

u 

velocity vector, m/s 

D h 

hydraulic diameter, m 

u 

velocity, m/s 

H 

height, m 

u 

average velocity, m/s 

L 

length, m 

A p 

pressure drop, Pa 

M 

amplitude of the mass flow rate, kg/s 

Pf 

Forchheimer coefficient, kg/m 4 

M g 

molar mass of the gas, kg/mol 

e 

porosity 

N 

number of pillars in the computational domain 

V 

regenerator effectiveness 

Nu Dh 

Nusselt number based on the hydraulic diameter 

K 

permeability, m 2 

0 

heat removal, J 

P 

dynamic viscosity, Pa s 

R 

ideal gas constant, J/(mol I<) 

9 

density, kg/m 3 

Re D„ 

Reynolds number based on the hydraulic diameter 

a 

electrical conductivity, S/m 

s 

Seebeck coefficient, V/K 

T 

viscous tensor, N/m 2 

S L 

horizontal pitch, m 



St 

vertical pitch, m 

Subscripts 

T 

absolute temperature, I< 

C 

cold side 

V 

volume, mm 3 

H 

hot side 

W 

work, J 

g 

gas 

^max 

vibration amplitude of the diaphragm, m 

loss 

loss 

ZT 

figure of merit 

s 

solid 

a fs 

solid surface area per unit volume, 1 /m 

out 

outlet 

c p 

specific heat, J/(kg K) 

r 

regenerator 

f 

frequency, Hz 



hfs 

convection coefficient, W/(m 2 K) 

Others 


k 

thermal conductivity, W/(m K) 

n 

unit vector normal to the boundary. 

m 

mass, kg 

V 

differential operator, 1 /m 

m 

mass flow rate, kg/s 




The traditional mechanical Stirling cooler uses two pistons to 
transfer the working fluid back and forth between hot and cold 
chambers that are separated by a regenerator. The heat is released 
from the hot chamber during the compression process and absor¬ 
bed from the cold chamber during the expansion process. Efforts on 
a prototype micro-scale Stirling cooler, which are documented in a 
series of cryocooler patents [12-15], addressed frictional losses and 
leakage concerns by replacing the conventional pistons and the 
associated linkages with electrostatically-driven diaphragms. 

Ceperley [16] recognized that sound waves can be used to 
replace the pistons or diaphragms for gas compression and 
displacement in the Stirling cooler. Thermoacoustic technology 
based on the Stirling cycle was thus developed. Reid et al. [17], Jin 
et al. [18], and Symko et al. [19] built thermoacoustic refrigeration 



Fig. 1. Theoretical COP from Eq. (1) as a function the temperature difference ( T H - T c ) 
for a thermoelectric cooler. 


systems and applied them in cryogenic cooling. Banjare et al. [20] 
and Zink et al. [21] performed CFD analysis to demonstrate the 
principle of the thermally-driven cooling. The temperature 
gradient achieved due to the thermoacoustic cooling, however, is 
limited by the critical temperature gradient [22]. 

In an earlier report, we presented the design concept of a new 
Stirling microcooler and a first-order thermodynamic (i.e., 
isothermal) analysis to evaluate the system performance [23]. A 
parametric study showed the effects of diaphragm phase lag, swept 
volume ratio between the hot space and cold space, and dead 
volume ratio on the cooling performance. For the Stirling micro¬ 
cooler, modeling challenges arise mainly from the complicated 
geometrical structures, (e.g., the large number of pillars in the 
regenerator) and complex dynamics (e.g., the motion of the di¬ 
aphragms). We herein address these issues. 

In this paper, more detailed numerical modeling that in¬ 
corporates compressible fluid flow, heat transfer, and solid me¬ 
chanics is used to study the performance of the full system. 
COMSOL [24], a multiphysics simulation software package, is used 
to perform the calculations. The rest of this paper is organized as 
follows. In Section 2, we present our design and a review of pre¬ 
vious work on modeling and applying the Stirling cycle, with a 
focus on sources of inefficiency and losses. To reduce the compu¬ 
tational time, in Section 3 the regenerator is isolated and its design 
is optimized. The full system performance is then evaluated in 
Section 4. 

2. Stirling cooler design, modeling, and inefficiencies 

We recently reported on the design of a new micro-scale Stirling 
cooler system, which includes two diaphragms and a regenerator 
that separates the hot and cold chambers [23]. By the design 
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requirements, each cooling element is 5 mm-long, 2.5 mm-wide, 
has a thickness of 0.15 mm, and is fabricated on a silicon wafer, as 
illustrated in Fig. 2(a) for one-half of an element. The components 
of one complete element are shown conceptually in Fig. 2(b). The 
design minimizes conduction heat losses across the 0.5 mm-long 
regenerator by distancing the compressor and expander assemblies 
with a low thermal conductivity passage, the regenerator, housing a 
working fluid (e.g., air) pressurized at 2 bar. The fluid flow is driven 
by electrostatically-actuated diaphragms. In the regenerator, an 
array of vertical silicon pillars serves as a thermal capacitor trans¬ 
ferring heat to and from the working gas during the cycle. Silicon 
pillars are also put in the space between the chambers and the 
regenerator (i.e., the dead space), for improving the heat transfer 
between the gas and the silicon around the chambers. The pillars in 
the dead space are part of one continuous piece of silicon, as 
illustrated in Fig. 2(a). Under operating conditions, the hot and cold 
diaphragms oscillate sinusoidally and 90° out of phase such that 
the heat is extracted to the cold chamber and released from the hot 
chamber. 

Classical theories, such as an isothermal model or an adiabatic 
model, can be applied to analyze Stirling coolers [25]. The 
isothermal analysis is based on the Schmidt cycle with sinusoidal 
volume variations. Assumptions are required to obtain a closed 
form-analytical expression for the cooling power. Notably, the gas 
temperature in the hot (cold) side is assumed to be the constant 
and to take on the temperature of the heat sink (source) and the 



(b) 


Fig. 2. Solid-model view of the Stirling microcooler element, (a) A single element is 
5 mm-long, 2.5 mm-wide, has a thickness of 150 gm, and is fabricated on a silicon 
wafer, (b) The assembled structure has five parts: the diaphragm layer in the middle, 
the top and bottom chamber substrates, and two sealing PDMS layers. 


regenerator is assumed to be perfect. The adiabatic analysis is more 
realistic as the working spaces tend to be adiabatic rather than 
isothermal in real machines. The adiabatic analysis assumes that 
the hot and cold spaces are thermally-insulated and that all the 
heat input and output occur at the heat exchangers. The gas flows 
out of the heat exchangers at the heat sink (source) temperature 
and then comes into the hot (cold) space. The regenerator is also 
assumed to be perfect. 

Martini [26], Wale et al. [27], Lee et al. [28], and Shoureshi [29] 
included more realistic losses and inefficiencies in the isothermal 
or adiabatic models. The important losses and inefficiencies in the 
Stirling cooler devices include: 1) Parasitic regenerator heat con¬ 
duction, whereby a heat flow develops along the device walls due 
to the temperature gradient between the hot and cold sides. 2) Due 
to the finite convection coefficient and limited heat capacity of the 
regenerator, when the gas flows from the hot side to the cold side, 
the output gas temperature is higher than the cold space temper¬ 
ature, thus reducing the maximum possible heat transfer. 3) The 
pressure drop when the gas flows through the regenerator. 4) The 
finite convection coefficients in the hot and cold exchangers. For 
different types of Stirling devices, losses and inefficiencies may also 
exist due to gas leakage, mechanical damping, and electrical effects 
[ 11 ]. 

In this study, the losses associated with the regenerator 
(parasitic heat flow, pressure drop, and insufficient heat transfer) 
and the finite heat transfer in the chambers are incorporated in a 
system-level model to evaluate the cooling performance of the 
device. Future modeling work may focus on the mechanical loss 
associated with driving the diaphragm and electrical losses from 
the circuits. 


3. Regenerator analysis 

Estimating the COP requires knowledge of the real work input 
and the heat transfer inefficiencies. In our design, inefficiencies 
come from the convective heat transfer resistance between the 
silicon substrate and the gas in the chamber, the insufficient heat 
transfer in the regenerator (i.e., the regenerator effectiveness is not 
unity), and the extra work required to overcome the pressure drop 
through the regenerator. The regenerator is thus a critical compo¬ 
nent of the Stirling cooler. Before performing the system-level 
modeling, we first model the regenerator to determine how its 
complicated geometry can be optimized. In this part of the analysis, 
the gas temperatures in the chambers are assumed to be the same 
as the heat source/sink temperatures. 

The purpose of the regenerator is to store and release heat from/ 
to the gas during the cycling. The regenerator design has the 
following requirements [11]: 1) A maximum ratio of the regener¬ 
ator heat capacity to the gas heat capacity. 2) A maximum heat 
transfer between the gas and the regenerator, requiring a large 
contact area. 3) A minimum pressure drop across the regenerator. 
4) Complete penetration of the heat in the regenerator material 
when it is heated or cooled. This last requirement can be achieved 
by using a solid material with small characteristic dimensions. As 
discussed in Ref. [23], circular silicon pillars with a diameter of 
20 pm are used for the regenerator solid structure due to the lim¬ 
itations of the heat penetration and available fabrication processes. 

Finite element analysis of two-dimensional fluid flow is used to 
analyze the arrangement of the pillars in the regenerator (of 
particular interest is the porosity) in order to balance the effects of 
the pressure drop and heat transfer. If inefficiencies reduce the heat 
transfer by Qioss and the work loss due to the pressure drop is Wi 0 S s, 
combining the ideal system cooling power and the work calculated 
by the isothermal model leads to a system COP of 
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hot-side surface 
< 



cold-side surface 

> 


n 


symmetric air silicon pillar 


Fig. 3. Two-dimensional regenerator computational domain and boundary conditions 
for the porosity optimization study (see also Table 1). The horizontal pitch is S L and the 
vertical pitch is S T . The regenerator length is L r . 


COP 


Q Qioss 

W + w loss • 



The computational domain and boundary conditions are illus¬ 
trated in Fig. 3. Taking advantage of symmetry, we only need to 
consider half of an array. The horizontal and vertical pitches are S L 
and St. The total length of the regenerator, L r , is 0.5 mm. The width 
of the studied section is SjI 2, which corresponds to a slice of the 
much wider regenerator. The air flows back and forth between the 
hot and the cold sides. The mass flow rate at the hot-side surface is 
taken to be m = Msin(2Tu/t), where M = 7zpfS T V H /2H r is the mass 
flow rate amplitude, H r is the regenerator height, and Vh is the 
volume of the hot gas space, which is 0.48 mm 3 . The volume of the 
chamber space is larger than the regenerator volume, which is 
0.17 mm 3 when the porosity of the regenerator is unity. The no-slip 
condition is applied at the gas-solid interfaces. The hot-side sur¬ 
face temperature condition is T = Th = 313.15 K if u n < 0 and 
VTn = 0 if u n > 0. The cold-side temperature condition is 
T= Tc = 288.15 K if u n < 0 and VT n = 0 if u n > 0. Here, u is the 
velocity vector and n is the unit vector normal to the boundary, as 
defined in Fig. 3. When the fluid flows out from the cold/hot side, 
the temperature is the heat source/sink temperature. Different 
porosities are obtained by adjusting the number of pillars, N, and 
their pitch. To assess the sensitivity of the results to the geometry, 
two SlISt ratios are considered, as shown in Table 1. e is the 
regenerator porosity. The time-averaged pressure drop across the 
regenerator per cycle, Ap , and the actual cold-side surface tem¬ 
perature, Tout, are obtained from the simulation results. The effec¬ 
tiveness of the regenerator is calculated from 


(Th — T ou t) 
(Th - T c ) ' 



3.1. Fluid compressibility 


Due to the pressure variations in the Stirling cooler (1.5— 
2.8 bar), fluid compressibility effects may be important. Most of the 
air properties relevant to this study do not change significantly 
within this pressure range (e.g., specific heat, thermal conductivity, 


Table 1 

Regenerator pillar configurations used in porosity optimization study. 


SlISt 

N 

S t /2 (pm) 

e 

0.289 

13 

65.7 

0.938 


15 

56.3 

0.916 


17 

49.3 

0.892 


19 

43.8 

0.864 


21 

39.4 

0.833 


23 

35.8 

0.798 

0.866 

8 

37.5 

0.933 


9 

32.8 

0.914 


10 

29.2 

0.892 


11 

26.3 

0.868 


12 

23.9 

0.842 


13 

21.9 

0.813 


and dynamic viscosity [30]). Only the density change is significant. 
Before performing the regenerator optimization, we first determine 
if compressibility effects need to be included. For the compressible 
flow, the air is assumed to be ideal, and the equation of state of a 
classical ideal gas is used to calculate density. The gas pressure as a 
function of time is calculated based on the isothermal model and is 
in the range of 1.5-2.8 bar. For the incompressible flow, the gas 
pressure is set to 2 bar. The gas pressure is specified at the cold-side 
surface of the computation domain. The operating-frequency- 
dependent pressure drop through the regenerator and the regen¬ 
erator effectiveness are plotted in Fig. 4(a) and (b) for both the 
compressible and incompressible flows. From the results, the effect 
of compressibility on the flow and heat transfer in the regenerator 
is small. We therefore use an incompressible flow computation for 
the optimization analysis of the regenerator. 

3.2. Porosity optimization 

With T h = 313.15 K and T c = 288.15 K, the Carnot COP is 11.5. The 
regenerator pressure drop and effectiveness as a function of the 
porosity and operating frequency are plotted in Fig. 5(a) and (b). 
The work loss per cycle is 

W loss = 2A p-V c . (4) 

The insufficient heat transfer loss per cycle is 

Qioss = (1 - V) x ( c p m) g X (Th - T c ). (5) 

Here, m is the maximum gas mass in the hot chamber. Using Eq. 
(2), the actual COP is plotted in Fig. 6(a) and (b) for two different 
pitch ratios. The results in Fig. 5 indicate that the pressure drop 
increases with decreasing porosity and increasing frequency, and 




f (Hz) 


Fig. 4. Frequency-dependence of (a) time-averaged pressure drop across the regen¬ 
erator and (b) regenerator effectiveness for SJS T = 0.289 and e = 0.892 and 0.833. The 
differences between considering compressible or incompressible flow are small. 
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(b) 



£ 


Fig. 5. Effect of operating frequency and porosity for Sl/St = 0.289 on (a) time- 
averaged pressure drop across the regenerator and (b) regenerator effectiveness. 


(a) 




Fig. 6. Predicted COP as a function of operating frequency and porosity for (a) SJ 
S T = 0.289 and (b) S L /S T = 0.866. 


that the heat transfer between the solid and the gas decreases with 
increasing porosity and frequency. An optimal porosity, as seen in 
Fig. 6, exists due to the large work loss at small porosities and the 
large heat transfer loss at high porosities. Even though the ar¬ 
rangements of the regenerator are significantly different (S^/ 
S T = 0.289 and 0.866), the optimal porosity is always near 0.9. 
Therefore, the key parameter in the regenerator that affects the 
system COP is the porosity. We note that the COP at Sl/St = 0.289 is 
slightly bigger than that at Si/Sj = 0.866 at higher frequency. As the 
gas temperatures in the chambers are assumed to be the same as 
the heat source/sink temperatures, meaning that the convection 
coefficients in the chambers are taken to be infinite, the COPs 
obtained here are overestimated. 


4. System evaluation 

In the previous section, the micro-Stirling cooler regenerator 
was modeled and we found that the optimal porosity is near 0.9. In 
this section, we choose a porosity of 0.892 at Si/Sj = 0.289 as a case 
study and evaluate the full system performance. It is not compu¬ 
tationally feasible to build a model that includes the full details of 
the regenerator geometry (i.e., all the pillars). To overcome the 
computational complexity brought about by the fine pillar struc¬ 
ture in the regenerator and in the dead space, a porous medium 
model was used to replace the pillars, allowing for a full-system 
finite element calculation. 


4.2. Porous medium governing equations 


As the Stirling cooler works by gas expansion and compression 
processes, the main physics in the system are the compressible 
laminar flow and heat transfer. In the porous medium region, a 
non-equilibrium model is used [31 ]. The governing equations for 
the gas phase are 


dp V-Qou) 
dt e 




p 

£ 


au 

at 


, N u 

+ (U-V)- 
£ 



V-T 

£ 


(£ + fr|u|)u, 



= V-(kVT) + + h fs a fs (T s - T). 

( 8 ) 

Here p is the pressure, p is the dynamic viscosity, k is the 
permeability, fa is the Forchheimer coefficient, T s is the solid tem¬ 
perature, hf s is the convection coefficient between the solid and the 
fluid, and a/ s is the solid surface area per unit volume, and f is the 
viscous tensor, 


P c p 


a t 

£ - + (u-V)T 


Tl 2 

Vu+ (Vll) J J --p(V-ll). 


T = p 

Assuming the gas to be ideal, the state equation is 


P 


pM, 


g 


RT ’ 


( 9 ) 

( 10 ) 


where R is the ideal gas constant and M g is the molar mass. 
In the solid phase, the energy equation is 


(i 




ke-WT s -h fs a fs {T s -T), 


( 11 ) 
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where p s and c Ps are the density and specific heat of the solid. k e is 
the effective thermal conductivity tensor of the solid phase in the 
porous medium. 


4.2. Specifying the permeability and Forchheimer coefficient 


The permeability k and Forchheimer coefficient ftp can be ob¬ 
tained based on the one-dimensional Darcy-Forchheimer equa¬ 
tion, which is [32] 


dp 

Ox 


-u + fipU 2 . 

K 


( 12 ) 



The relationship between the pressure drop and the velocity for the 
regenerator with pillars is obtained by steady-state simulation and is 
plotted in Fig. 7. The permeability and Forchheimer coefficient are then 
calculated by the least squares method according to Eq. (12), which is 
also plotted in Fig. 7. The best-fit permeability is 9.77 x 1CT 11 m 2 
and the best-fit Forchheimer coefficient is 9823 kg/m 4 . 

4.3. Convection coefficient 

The convection coefficient of the regenerator is also evaluated 
from steady-state simulations. The Nusselt number and Reynolds 


number are defined as 


hfrDl ? 

Nu Dll = , 

(13) 

Dp f ,m h 

Re °‘ - n ’ 

(14) 


Fig. 8. Regenerator Nusselt number as a function of Reynolds number for the current 
calculations and from existing correlations [33-35]. The current simulation results are 
used for the system modeling. 


and the dead-space structure. For the solid phase of the porous 
medium in the dead space, an isotropic thermal conductivity 
(1 - c)k S iiicon is assumed. The thermal conductivity of the solid 
phase in the regenerator is anisotropic due to the vertical array of 
silicon pillars, as illustrated in Fig. 2(a). The out-of-plane thermal 
conductivity is the same as that of the porous medium of the 
dead space. The in-plane thermal conductivity is taken to be that 
of air. 

As the diaphragms are driven electrostatically, their actual 
motions are complicated and require detailed study. To simplify the 
model and the analysis, a sinusoidal motion with a 90° phase lag is 
applied for the cold and hot diaphragms. The displacement of the 
cold diaphragm is 


where Dh = 4e\af S is the hydraulic diameter and u is the average 
velocity through the regenerator. The flow in the system is laminar 
as the Reynolds number in the regenerator is smaller than 200. The 
Nusselt number is calculated using the log-mean temperature dif¬ 
ference when a constant temperature is applied to the pillar walls 
and is plotted in Fig. 8. Fitting the raw data gives Nu Dh = 13.389Re D/7 
which is used for the system modeling. Some available experi¬ 
mental correlations of the Nusselt number for a tube bank are 
plotted in Fig. 8 [33-35]. These correlations focus on the flow in a 
larger Reynolds number range (Re Dh > 300 in the present study). 
The convection coefficient is also influenced by the arrangement of 
pillars in the regenerator. Thus, differences between the calculation 
results and correlations are observed. 

4.4. Computational domain and boundary conditions 

The computational domain for the full system model is shown 
in Fig. 9. We use the porous medium to replace the regenerator 



Fig. 7. Pressure drop across the regenerator as a function of average fluid velocity 
between (i) the full pillar structure and (ii) the porous medium model. 


z(x,y,t) = Z maX ( X j,)Sin(27r/t). (15) 

The displacement of the hot diaphragm is 

z(x,y,t) = Z max(x , y) sin(27t/t-|), (16) 

where /is the operating frequency and Z ma x(x,y) is the vibration 
amplitude of the diaphragm. We assume that the gas space in the 
chamber is a spherical cap. The cap base is a circle with a 2.25 mm 
diameter and the maximum height of the cap is 120 pm. So, 


5 mm 





L . 

3.5 mm 


A’ 


T c cold diaphragm 


hot diaphragm 


H 



silicon 


A-A’ 0.15 mm 

porous medium in the dead space 


air 


porous medium in the regenerator 


Fig. 9. Computational domain and geometry for the three-dimensional full-system 
model. The hot and cold diaphragms, which have a 90° phase lag, prescribe the moving 
boundaries of the fluid flow. The regenerator and the dead space are replaced by the 
porous medium. 
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Fig. 10. Space-averaged cooling capacity as a function of time at an operating fre¬ 
quency of 100 Hz when the porosity is 0.892. Steady state is achieved after four cycles. 


Zmax(Xy) represents the cap height at different positions. The hot 
and cold diaphragms prescribe the moving boundaries of the 
fluid flow. The arbitrary Lagrangian-Eulerian (ALE) moving-mesh 
method was used to handle the gas domain where there is a mesh 
deformation. The no-slip condition is applied to the solid-gas 
interfaces. The initial pressure in the system is 2 bar. The tem¬ 
perature of the end surface of the cold side is constant and is 
Tc = 288.15 K. The end surface in the hot side has a constant 
temperature of Th = 313.15 K. The remaining surfaces are ther¬ 
mally insulated. The initial temperature for the whole region is 
293.15 K. The number of mesh elements for the system model is 
347,773. Parallel computation was performed by using 4 nodes, 
each of which has 16 CPUs and 16 GB of memory. Each 


computation took about 20 h. The solver used for the computa¬ 
tion is MUMPS [36]. 


4.5. Results and discussions 

The effect of the operating frequency on the cooling capacity 
and COP is now explored. The frequency range is 100-800 Hz. 
Steady-state is reached after four cycles of the diaphragm operation 
at an operating frequency of 100 Hz, as shown in Fig. 10. At steady- 
state, the temperature contours of the gas in the system and the 
silicon around the chamber at one quarter and three quarters of one 
cycle are shown in Fig. 11(a) and (b) when the frequency is 100 Hz. 
In Fig. 11(a), the cold diaphragm is completely compressed and the 
hot gas space has its highest temperature at this point in the cycle. 
During the compression process, the heat is released from the hot 
gas and transferred to the heat sink. In Fig. 11(b), the cold gas space 
is expanded. At this time, the cold gas space has its lowest tem¬ 
perature. This process is an expansion in which the heat from the 
heat source is absorbed by the gas in the cold space. Through cy¬ 
clings, the heat is transferred to the heat sink and absorbed from 
the heat source continually, as we require. The temperature con¬ 
tours at one quarter and three quarters of one cycle at a frequency 
of 800 Hz are shown in Fig. 11(c) and (d). The heat transfer process 
is similar compared to the results at 100 Hz. However, the 
maximum temperature in the hot gas space is much higher and the 
minimal temperature in the cold gas is much lower at 800 Hz. The 
average cold and hot gas space temperatures at different fre¬ 
quencies are plotted in Fig. 12. The cold gas temperature is colder 



Fig. 11. Temperature contours for the full Stirling cooler system simulation at: (a) One quarter of the cycle at an operating frequency of 100 Hz. (b) Three quarters of the cycle at an 
operating frequency of 100 Hz. (c) One quarter of the cycle at an operating frequency of 800 Hz. (d) Three quarters of the cycle at an operating frequency of 800 Hz. 
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Fig. 12. Time- and space-averaged gas temperature in the cold and hot chambers after 
the system comes to the steady state as a function of frequency. 



q (W/cm 2 ) 

Fig. 14. Multiphysics simulation prediction of the system COP as a function of the 
cooling capacity. The theoretical COP is 11.5 with T H = 313.15 K and T c = 288.15 K. 


when the frequency is higher, and the hot gas temperature has the 
opposite trend. This result indicates that the heat transfer resis¬ 
tance between the gas and the silicon around the chamber affects 
the system performance, and this heat transfer resistance increases 
with an increase of the frequency. 

The heat flux coming into the cold side from the end surface 
represents the cooling capacity of the Stirling cooler. The time- 
dependent pressure and volume variations in the element are ob¬ 
tained from this simulation. The PV work can therefore be calcu¬ 
lated and the COP can be evaluated. The results of cooling capacity 
are shown in Fig. 13. The theoretical cooling capacities of the 
isothermal model and the adiabatic model are also provided in 
Fig. 13. The adiabatic model result is closer to the simulation as the 
isothermal assumption in chambers is not good due to the high 
operating frequency, as shown in Fig. 11. The cooling capacity in¬ 
creases almost linearly from 0.7 to 5.1 W/cm 2 when the frequency 
increases from 100 Hz to 800 Hz. The COP with respect to the 
cooling capacity is plotted in Fig. 14. The COP decreases from 7.4 to 
2.4 when the operating frequency and the cooling capacity in¬ 
crease. The reason is that even though the cooling power increases 
when the frequency increases, the thermodynamic work in the 
system increases more quickly, making the COP lower. An appro¬ 
priate frequency should be considered according to the require¬ 
ment of the real situation. When the operating frequency is 600 Hz, 
the COP is 2.93, which is about 25% of the Carnot COP, and the 
cooling capacity is 4.2 W/cm 2 . In a standard design of a thermo¬ 
electric cooler (i.e., the material thickness is larger than 1 mm), the 
maximum cooling capacity is about 2 W/cm 2 and the best COP is 
2.85 when the ZT of the material is 1 and the temperature differ¬ 
ence is 15 I< [2]. 



Fig. 13. System cooling capacity as a function of the frequency from the multiphysics 
simulations, the isothermal model, and the adiabatic model. With an increase of the 
operating frequency, the cooling capacity increases. The adiabatic model result is closer 
to the simulation result than the isothermal model. 


5. Conclusions 

In this paper, a new Stirling micro-scale cooler has been 
modeled and simulated numerically. The computations of multi¬ 
physics processes, incorporating compressible fluid flow, heat 
transfer, porous medium, solid mechanics, and moving mesh, have 
been successfully implemented. 

The regenerator performance in the Stirling cooler was studied 
and the effects of pressure drop and heat transfer were analyzed. 
Optimizations indicated that the optimal porosity of the regener¬ 
ator is near 0.9. The complete system modeling predicted the 
system-level thermal performance. Parametric studies of the 
design demonstrated the effect of the operating frequency on the 
cooling capacity and the COP of the system. When the system is 
operated at 600 Hz, the cooling capacity is 4.2 W/cm 2 and the 
system COP is 2.93 when the temperatures of the heat sink and heat 
source are 313.15 I< and 288.15 K. 

Acknowledgments 

The authors express their appreciation to Dr. Suresh Santhanam 
for helpful discussions during the research. This work was sup¬ 
ported by the Defense Advanced Research Projects Agency (DARPA) 
and the U.S. Army Aviation and Missle Research, Development, and 
Engineering Center (AMRDEC) under Grant No. W31P4Q-10-1- 
0015. The views and conclusions contained in this document are 
those of the authors and should not be interpreted as representing 
the official policies, either expressed or implied, of DARPA, the U.S. 
Army, or the U.S. Government. 

References 

[1] M.E. Moran, D.M. WesolekO, B.T. Berhane, K.J. Rebello, Microsystem cooler 
development, in: 2nd International Energy Conversion Engineering Confer¬ 
ence, Providence, RI, 2004. 

[2] J. Sharp, J. Bierschenk, H.B. Lyon, Overview of solid-state thermoelectric re¬ 
frigerators and possible applications to on-chip thermal management, in: 
Proceedings of the IEEE vol. 94, 2006, pp. 1602—1612. 

[3] L.W. da Silva, M. Kaviany, Fabrication and measured performance of a first- 
generation microthermoelectric cooler, Journal of Microelectromechanical 
Systems 14 (5) (2005) 1110-1117. 

[4] R. Venkatasubramanian, E. Siivola, T. Colpitts, B. O’Quinn, Thin-film thermo¬ 
electric devices with high room-temperature figures of merit, Nature 413 
(2001) 597-602. 

[5] C. LaBounty, A. Shakouri, J.E. Bowers, Design and characterization of thin film 
microcoolers, Journal Applied Physics 89 (7) (2001) 4059—4064. 

[6] A. Majumdar, Thermoelectricity in semiconductor nanostructures, Science 
303 (5659) (2004) 777-778. 

[7] W. Kim, J. Zide, A. Gossard, D. Klenov, S. Stemmer, A. Shakouri, A. Majumdar, 
Thermal conductivity reduction and thermoelectric figure of merit increase by 
embedding nanoparticles in crystalline semiconductors, Physical Review 
Letters 96 (2006) 045901. 
































52 


D. Guo et al. / International Journal of Thermal Sciences 74 (2013) 44—52 


[8] Y.S. Ju, Solid-state refrigeration based on the electrocaloric effect for elec¬ 
tronics cooling, Journal of Electronic Packaging 132 (2010) 041004. 

[9] M. Valant, Electrocaloric materials for future solid-state refrigeration tech¬ 
nologies, Progress in Materials Sciences 57 (6) (2012) 980—1009. 

[10] S.L. Garrett, Reinventing the engine, Nature 339 (1999) 303—305. 

[11] J.F. Burger, Cryogenic Microcooling. PhD thesis, University ofTwente, 2001. 

[12] N.B. Stetson, Miniature integral Stirling cryocooler, U.S. Patent No. 5,056,317, 
1991. 

[13] R. Solomon, Integrated refrigerated computer module, U.S. Patent No. 
5,349,823, 1994. 

[14] L. Bowman, D.M. Berchowitz, I. Uriell, Microminiature Stirling cycle cry- 
ocoolers and engines, U.S. Patent No. 5,457,956, 1995. 

[15] M.E. Moran, Micro-scalable thermal control device, US Patent No. 6,385,973 
Bl, 2002. 

[16] P.H. Ceperley, Gain and efficiency of a short traveling wave heat engine, 
Journal of Acoustical Society America 77 (3) (1985) 1239—1244. 

[17] R.S. Reid, W. Ward, G. Swift, Cyclic thermodynamics with open flow, Physical 
Review Letters 80 (21) (1998) 4617-4620. 

[18] T. Jin, G.B. Chen, B.R. Wang, S.Y. Zhang, Application of thermoacoustic effect to 
refrigeration, Review of Scientific Instruments 74 (1) (2003) 677—679. 

[19] O.G. Symko, E. Abdel-Rahman, Y.S. Kwon, M. Emmi, R. Behunin, Design and 
development of high-frequency thermoacoustic engines for thermal man¬ 
agement in microelectronic, Microelectronic Journal 35 (2) (2004) 185—191. 

[20] Y.P. Banjare, R.K. Sahoo, S.K. Sarangi, CFD simulation inertance tube pulse tube 
refrigerator, in: 19th National and 8th ISHMT-ASME Heat and Mass Transfer 
Conference, 2008. Paper # EXM-7. 

[21] F. Zink, J. Vipperman, L. Schaefer, CFD simulation of thermoacoustic 
cooling, International Journal of Heat and Mass Transfer 53 (19) (2010) 
3940-3946. 

[22] G.W. Swift, Thermoacoustic engines, Journal of Acoustical Society America 84 
(4) (1988) 1145-1180. 

[23] D. Guo, J. Gao, A.J.H. McGaughey, G.K. Fedder, M.E. Moran, S.C. Yao, Design and 
evaluation of a MEMS-based Stirling microcooler, Journal of Heat Transfer, in 
press. 


[24] D. Guo, A.J.H. McGaughey, G.K. Fedder, M. Lee, S.C. Yao, Modeling system 
dynamics in a MEMS-based Stirling cooler, in: COMSOL Conference, Boston, 
MA, 2011. 

[25] G. Walker, Stirling Cycle Machines, Clarendon Press, Oxford, 1976. 

[26] W.R. Martini, Thermodynamic design of Stirling engines by computer, Martini 
Engineering, Richland, Washington, 1980. 

[27] E.R. Wale, J.I. Smith Jr., A mathematical model for steady operation of Stirling- 
type engines, Journal of Engineering for Gas Turbines and Power 90 (1) (1968) 
45-50. 

[28] K. Lee, I.P. Krepchin, W.M. Toscano, Thermodynamic description of an adia¬ 
batic second order analysis for Stirling engines, in: Proceedings of the 16th 
IECEC, 1981, pp. 1919-1924. Paper # 819794. 

[29] R. Shoureshi, Analysis and Design of Stirling Engines for Waste-Heat Recovery. 
PhD thesis, MIT, 1981. 

[30] K. Kadoya, N. Matsunaga, A. Nagashima, Viscosity and thermal conductivity of 
dry air in the gaseous phase, Journal of Physical and Chemical Reference Data 
14 (4) (1985) 947-970. 

[31] R.C. Tew, T. Simon, D. Gedeon, M. Ibrahim, W. Rong, An initial non¬ 
equilibrium porous media model for CFD simulation of Stirling regenerators, 
in: 4th International Energy Conversion and Engineering Conference, San 
Diego, CA, 2006. 

[32] A. Amiri, K. Vafai, Transient analysis of incompressible flow through a packed 
bed, International Journal of Heat and Mass Transfer 41 (24) (1998) 4259— 
4279. 

[33] A.A. Zukauskas, Heat transfer from tubes in cross flow, in: Advances in Heat 
Transfer, Academic Press, New York, 1972, pp. 93—160. 

[34] S. Whitaker, Forced convection heat-transfer correlations for flow in pipes, 
past flat plates, single cylinders, single spheres, and for flow in packed-beds 
and tube bundles, AIChE Journal 18 (2) (1972) 361—371. 

[35] T.H. Hwang, S.C. Yao, Crossflow heat transfer in tube bundles at low Reynolds 
numbers, Journal of Heat Transfer 108 (3) (1986) 697—700. 

[36] P.R. Amestoy, I.S. Duff, J.-Y. L’Excellent, Multifrontal parallel distributed 
symmetric and unsymmetric solvers, Computer Methods in Applied Me¬ 
chanics and Engineering 184 (2000) 501—520. 


